home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / sign_test.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  209 lines
  1. ; $Id: sign_test.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7. function gaussint1,x
  8. ; gaussinit1 returns the probabilty of obtaining x or something
  9. ; more extreme.
  10.  
  11.   if x le 0 then return,gaussint(x)          $
  12.  else return,1 - gaussint(X)
  13.  
  14. end 
  15.  
  16.  
  17. pro sign_test,Data, Diff,Prob,Names=Names,List_Name=Ln,   $
  18.               Missing=M,NoPrint=NP
  19. ;+ 
  20. ; NAME:
  21. ;    SIGN_TEST
  22. ;
  23. ; PURPOSE:
  24. ;    To test the null hypothesis that two populations have the same 
  25. ;    distribution -i.e. F(x) = G(x) against the alternative that their 
  26. ;    distributions differ in location- i.e F(x) = G(x+a).Sign_test 
  27. ;    pairwise tests the populations in Data.
  28. ;  
  29. ; CATEGORY:
  30. ;    Statistics.
  31. ;
  32. ; CALLING SEQUENCE:
  33. ;    SIGN_TEST, Data, [Diff,Prob,Names=Names]
  34. ;
  35. ; INPUTS:
  36. ;    Data:    Two-dimensional array. Data(i,j) = the jth observation from
  37. ;        the ith population.
  38. ;
  39. ; KEYWORDS:  
  40. ;    NAMES:    Vector of user supplied names for the populations to be 
  41. ;        used in the output.
  42. ;
  43. ;    LIST_NAME:    Name of output file. Default is to the screen.
  44. ;
  45. ;      MISSING:    Value used as a place holder for missing data.  Pairwise 
  46. ;        handling of missing data.
  47. ;           
  48. ;      NOPRINT:    Flag, if set, to suppress output to the screen.
  49. ;
  50. ; OUTPUT:
  51. ;    Table written to the screen showing for each pair of populations 
  52. ;    the number of positive differences in observations.  Also, table of 
  53. ;     probabilties for each population pair giving the two-tailed 
  54. ;    significance of the results in the first table.
  55. ;
  56. ;OPTIONAL OUTPUT PARAMETERS: 
  57. ;    Diff:    Two-dimensional array of positive differences.
  58. ;        Diff(i,j) = number of observations in population
  59. ;        i greater than the corresponding observation in population j.
  60. ;
  61. ;    Prob:    Two-dimensional array. Prob(i,j) = probability of 
  62. ;        Diff(i,j) or something more extreme.
  63. ;                          
  64. ;RESTRICTIONS:
  65. ;      All populations have the same sample size.
  66. ;
  67. ;COMMON BLOCKS: 
  68. ;     None.
  69. ;
  70. ;PROCEDURE:
  71. ;    For each pair of populations, the diffence between corresponding
  72. ;    observations is computed and a count is made of the positive and
  73. ;    negative differences.  The probability of the count is computed
  74. ;    under the assumption that the distributions are the same - i.e.
  75. ;    the probability of a negative difference = the probability of a 
  76. ;    positive difference = .5.  For sample size > 25, the binomial 
  77. ;    distribution is approximated with a normal distribution for computing
  78. ;    Prob.
  79. ;-
  80.  
  81.  
  82. On_Error,2
  83. SD= size(Data)
  84.  
  85. if( N_Elements( Ln) NE 0) THEN openw,unit,/get,Ln else unit=-1
  86.  
  87. if ( SD(0) NE 2) THEN BEGIN
  88.    printf,unit, 'sign_test- Data array has wrong dimension'
  89.    goto, DONE
  90. ENDIF
  91.  
  92.  
  93. C=SD(1)
  94. R= SD(2)
  95.  
  96.  
  97.  
  98. Diff = Fltarr(C,C) 
  99. Prob = Diff +1.0
  100.  
  101.  for i = 0l,C-2 DO  BEGIN
  102.  
  103.    D1 = Replicate(1.0,C-i-1) # Data(i,*)  ;compute differences
  104.    Temp = Data(i+1:*,*) - D1
  105.  
  106.    if (N_Elements(M) NE 0) THEN BEGIN   ; Handle missing data
  107.      here = where( Data(i+1:*,*) EQ M, count) 
  108.    if count  NE 0 THEN Temp(here)=0
  109.    here = where(D1 EQ M,count)
  110.    if count NE 0 THEN  Temp( here) = 0
  111.    ENDIF
  112.  
  113.    here = where (Temp NE 0,count)
  114.    if ( count NE 0) THEN BEGIN    
  115.      Temp1 = Temp
  116.      Temp1 ( here) = 1   ; If diff =0 then discard observation
  117.      PopSize = Temp1 # Replicate(1,R)   
  118.                   ; compute number of observation per column
  119.  
  120.  
  121.  
  122.      here = where(Temp LE 0, count)     ;count positives
  123.      if (count ne 0) THEN  $
  124.         Temp(here) = 0                 
  125.      here = where(Temp NE 0,count)
  126.      if ( count NE 0) THEN $
  127.      Temp(here) = 1
  128.      PosNo = Temp # Replicate(1,R)     
  129.  
  130.  
  131.  
  132.      Diff(i+1:*,i) = PosNo
  133.  
  134.  
  135.      for j =long(i+1),C-1 DO BEGIN
  136.          k=j-i-1
  137.         if Popsize(k) eq 0 THEN BEGIN
  138.            printf,unit,"sign_test- Data are all the same or missing"
  139.            printf,unit,"           for columns ",i, " and ",j
  140.            Diff(i,j) = -1 & Diff(j,i) = -1
  141.            Prob(i,j) = -1 & Prob(j,i) =-1
  142.         ENDIF ELSE $
  143.         if PopSize(k) GT 25 THEN      $       
  144.                           ;approximate binomial with normal
  145.            Prob(j,i) =2*      $
  146.                  Gaussint1((2*Diff(j,i)-       $
  147.                     PopSize(k))/sqrt(PopSize(k)))  $
  148.  
  149.          else if Diff(j,i) GT PopSize(k)/2 THEN   $          
  150.                  prob(j,i) =      $
  151.                     2*binomial(Diff(j,i),PopSize(k) ,.5)     $
  152.          else if Diff(j,i) eq PopSize(k)/2 THEN       $ 
  153.              prob(j,i) = 1           $
  154.              else  prob(j,i) =        $
  155.                    2*(1- binomial(Diff(j,i)+1,PopSize(k),.5))
  156.    
  157.       ENDFOR
  158.  
  159.  
  160.      Diff(i,i+1:*) = PopSize - Diff(i+1:*,i)
  161.      Prob(i,i+1:*) = Prob(i+1:*,i)
  162.  
  163.  
  164.  ENDIF ELSE BEGIN
  165.         printf,unit,'sign_test- all data are missing for column ',i
  166.         printf,unit,'           or data the same in all columns after column',i
  167.         printf,unit," "
  168.         Diff(i+1:*,i) = -1 & Diff(i,i+1:*) = -1
  169.         Prob(i,i+1:*) = -1 & Prob(i+1:*,i) = -1
  170.         ENDELSE
  171. ENDFOR
  172.  
  173.  SN =Size(Names)
  174.  if (SN(1) EQ 0) THEN BEGIN
  175.      I = INDGEN(C)
  176.      Names=['pop'+StrTrim(I,2)]  
  177.  ENDIF ELSE                  $
  178.    if ( SN(1) LT C) THEN BEGIN
  179.      I = Indgen(C)
  180.      printf,unit,'sign_test- missing names'
  181.      Names=[Names, 'pop'+StrTrim(I(SN(1):C-1),2)]
  182.   ENDIF
  183.  
  184.  
  185.  
  186. if( Not keyword_set(NP)) THEN BEGIN
  187.  
  188.   printf,unit, " Table of Count Differences"
  189.   printf,unit, " "
  190.   printf,unit, format ='(8X,16(A8,2x))', Names
  191.  
  192.   for i= 0,C-1 do                       $
  193.      printf,unit, format='(A8,16(I8,2X))',Names(i),Diff(*,i)
  194.  printf,unit, " "
  195.  printf,unit, "Table of Probabilities:"
  196.  printf,unit," "
  197.  printf,unit, format ='(8X,16(A10,2x))', Names
  198.  
  199.  for i= 0,C-1 do                      $
  200.      printf,unit, format='(A8,16(G10.5,2X))',Names(i),  $
  201.                 Prob(*,i)
  202.  ENDIF
  203.  
  204.  
  205. DONE:
  206.    if ( unit NE -1) THEN Free_Lun,unit
  207.    RETURN
  208.    END
  209.